plus sc10x
run cellranger with Spp1-EGFP added reference index
load dependancies
Spp1-GFP mice, SIMPLE induced stroke,
P28 7d: GFP(Spp1+) x2, MIG(Microglia) x2, CTR(Contra) x1
P05 3d: GFP(Spp1+) x2, MIG(Microglia) x2, CTR(Contra) x1
loading 70k cells, CX3CR1+
demo, cellranger called 44,614 cells
plus with Spp1-EGFP, cellranger called 47,886 cells
load 10x data
filt.10x <- Read10X(data.dir = "../output_plus_exogene/filtered_feature_bc_matrix/")## 10X data contains more than one type and is being returned as a list containing matrices of each type.GEX <- filt.10x$`Gene Expression`
FB <- filt.10x$`Antibody Capture`check datasets
dim(GEX)## [1] 32286 47886GEX[1:6,1:6]## 6 x 6 sparse Matrix of class "dgCMatrix"
## AAACCCAAGAAATCCA-1 AAACCCAAGAATCGTA-1 AAACCCAAGAGCCTGA-1
## Xkr4 . . .
## Gm1992 . . .
## Gm19938 . . .
## Gm37381 . . .
## Rp1 . . .
## Sox17 . . .
## AAACCCAAGATGCTAA-1 AAACCCAAGCACCTGC-1 AAACCCAAGCGATCGA-1
## Xkr4 . . .
## Gm1992 . . .
## Gm19938 . . .
## Gm37381 . . .
## Rp1 . . .
## Sox17 . . .dim(FB)## [1] 10 47886FB[,1:6]## 10 x 6 sparse Matrix of class "dgCMatrix"
## AAACCCAAGAAATCCA-1 AAACCCAAGAATCGTA-1 AAACCCAAGAGCCTGA-1
## P28.GFP.1 12 49 392
## P28.GFP.2 26 10 29
## P28.MIG.1 11 9 21
## P28.MIG.2 5 1 1
## P28.CTR.1 4 . 1
## P05.GFP.1 139 . 4
## P05.GFP.2 15 3 7
## P05.MIG.1 2 4 2
## P05.MIG.2 2 2 209
## P05.CTR.1 5 2 4
## AAACCCAAGATGCTAA-1 AAACCCAAGCACCTGC-1 AAACCCAAGCGATCGA-1
## P28.GFP.1 3 23 10
## P28.GFP.2 28 68 6
## P28.MIG.1 6 23 9
## P28.MIG.2 51 7 9
## P28.CTR.1 . 5 .
## P05.GFP.1 4 10 3
## P05.GFP.2 3 266 4
## P05.MIG.1 3 78 4
## P05.MIG.2 4 259 188
## P05.CTR.1 40 8 4rowSums(FB)## P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1 P05.GFP.2 P05.MIG.1
## 1648415 4014042 1456711 974485 312662 1004996 1222480 876714
## P05.MIG.2 P05.CTR.1
## 1141690 824693rowSums(FB>0)## P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1 P05.GFP.2 P05.MIG.1
## 47701 47859 47690 46474 35812 46092 47378 44721
## P05.MIG.2 P05.CTR.1
## 46311 42039FB
# build seurat object
FB.seur <- CreateSeuratObject(counts = FB,project = "plus_SIMPLE")
FB.seur## An object of class Seurat
## 10 features across 47886 samples within 1 assay
## Active assay: RNA (10 features, 0 variable features)rownames(FB.seur)## [1] "P28.GFP.1" "P28.GFP.2" "P28.MIG.1" "P28.MIG.2" "P28.CTR.1" "P05.GFP.1"
## [7] "P05.GFP.2" "P05.MIG.1" "P05.MIG.2" "P05.CTR.1"normalization
perform Seurat ’Demultiplexing with hashtag oligos(HTOs)
ref https://satijalab.org/seurat/articles/hashing_vignette.html
# normalize
FB.seur <- NormalizeData(FB.seur)
FB.seur <- FindVariableFeatures(FB.seur, selection.method = "mean.var.plot")
FB.seur <- ScaleData(FB.seur, features = VariableFeatures(FB.seur))## Centering and scaling data matrix# independent assay for HTO data
FB.seur[['HTO']] <- CreateAssayObject(counts = FB)
FB.seur <- NormalizeData(FB.seur, assay = 'HTO', normalization.method = 'CLR')## Normalizing across featurescheck demux
first define FB colors based on conditions
scales::show_col(ggsci::pal_igv("default")(49))
color.FB <- ggsci::pal_igv("default")(49)[c(2,13,33,1,15,
34,26,28,40,18)]
level.FB <- c(rownames(FB.seur))
color.FBraw <- c("#FF6C91","lightgrey",color.FB)scales::show_col(color.FBraw, ncol = 5)
scales::show_col(color.FB, ncol = 5)
par(mfrow=c(2,5))
#plot(sort(rowSums(t(FB))),
# xlab="reordered cells",
# ylab="UMI counts", log="xy",
# main="sum")
plot(sort(t(FB)[,1]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[1])
plot(sort(t(FB)[,2]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[2])
plot(sort(t(FB)[,3]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[3])
plot(sort(t(FB)[,4]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[4])
plot(sort(t(FB)[,5]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[5])
plot(sort(t(FB)[,6]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[6])
plot(sort(t(FB)[,7]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[7])
plot(sort(t(FB)[,8]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[8])
plot(sort(t(FB)[,9]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[9])
plot(sort(t(FB)[,10]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[10])
range ‘q’
q0.50 ~ 0.99
FB.test <- FB.seur
#cut.list <- list()
table.demux <- data.frame(sample=c("quantile","Doublet",rownames(FB.test),"Negative"))
for(i in 50:99){
FB.test <- HTODemux(FB.test, assay = "HTO", positive.quantile = i/100)
table.demux <- cbind(table.demux,
c(i/100,table(FB.test$hash.ID)[c("Doublet",rownames(FB.test),"Negative")]))
}## Cutoff for P28.GFP.1 : 7 reads## Cutoff for P28.GFP.2 : 15 reads## Cutoff for P28.MIG.1 : 7 reads## Cutoff for P28.MIG.2 : 4 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 3 reads## Cutoff for P05.GFP.2 : 7 reads## Cutoff for P05.MIG.1 : 2 reads## Cutoff for P05.MIG.2 : 4 reads## Cutoff for P05.CTR.1 : 2 reads## Cutoff for P28.GFP.1 : 7 reads## Cutoff for P28.GFP.2 : 15 reads## Cutoff for P28.MIG.1 : 8 reads## Cutoff for P28.MIG.2 : 4 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 3 reads## Cutoff for P05.GFP.2 : 7 reads## Cutoff for P05.MIG.1 : 2 reads## Cutoff for P05.MIG.2 : 5 reads## Cutoff for P05.CTR.1 : 2 reads## Cutoff for P28.GFP.1 : 7 reads## Cutoff for P28.GFP.2 : 16 reads## Cutoff for P28.MIG.1 : 8 reads## Cutoff for P28.MIG.2 : 4 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 3 reads## Cutoff for P05.GFP.2 : 7 reads## Cutoff for P05.MIG.1 : 2 reads## Cutoff for P05.MIG.2 : 5 reads## Cutoff for P05.CTR.1 : 2 reads## Cutoff for P28.GFP.1 : 7 reads## Cutoff for P28.GFP.2 : 16 reads## Cutoff for P28.MIG.1 : 8 reads## Cutoff for P28.MIG.2 : 4 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 3 reads## Cutoff for P05.GFP.2 : 7 reads## Cutoff for P05.MIG.1 : 2 reads## Cutoff for P05.MIG.2 : 5 reads## Cutoff for P05.CTR.1 : 2 reads## Cutoff for P28.GFP.1 : 7 reads## Cutoff for P28.GFP.2 : 16 reads## Cutoff for P28.MIG.1 : 8 reads## Cutoff for P28.MIG.2 : 4 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 3 reads## Cutoff for P05.GFP.2 : 7 reads## Cutoff for P05.MIG.1 : 2 reads## Cutoff for P05.MIG.2 : 5 reads## Cutoff for P05.CTR.1 : 3 reads## Cutoff for P28.GFP.1 : 8 reads## Cutoff for P28.GFP.2 : 16 reads## Cutoff for P28.MIG.1 : 8 reads## Cutoff for P28.MIG.2 : 4 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 3 reads## Cutoff for P05.GFP.2 : 8 reads## Cutoff for P05.MIG.1 : 3 reads## Cutoff for P05.MIG.2 : 5 reads## Cutoff for P05.CTR.1 : 3 reads## Cutoff for P28.GFP.1 : 8 reads## Cutoff for P28.GFP.2 : 16 reads## Cutoff for P28.MIG.1 : 9 reads## Cutoff for P28.MIG.2 : 4 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 3 reads## Cutoff for P05.GFP.2 : 8 reads## Cutoff for P05.MIG.1 : 3 reads## Cutoff for P05.MIG.2 : 5 reads## Cutoff for P05.CTR.1 : 3 reads## Cutoff for P28.GFP.1 : 8 reads## Cutoff for P28.GFP.2 : 17 reads## Cutoff for P28.MIG.1 : 9 reads## Cutoff for P28.MIG.2 : 5 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 3 reads## Cutoff for P05.GFP.2 : 8 reads## Cutoff for P05.MIG.1 : 3 reads## Cutoff for P05.MIG.2 : 5 reads## Cutoff for P05.CTR.1 : 3 reads## Cutoff for P28.GFP.1 : 8 reads## Cutoff for P28.GFP.2 : 17 reads## Cutoff for P28.MIG.1 : 9 reads## Cutoff for P28.MIG.2 : 5 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 3 reads## Cutoff for P05.GFP.2 : 8 reads## Cutoff for P05.MIG.1 : 3 reads## Cutoff for P05.MIG.2 : 6 reads## Cutoff for P05.CTR.1 : 3 reads## Cutoff for P28.GFP.1 : 8 reads## Cutoff for P28.GFP.2 : 17 reads## Cutoff for P28.MIG.1 : 9 reads## Cutoff for P28.MIG.2 : 5 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 4 reads## Cutoff for P05.GFP.2 : 8 reads## Cutoff for P05.MIG.1 : 3 reads## Cutoff for P05.MIG.2 : 6 reads## Cutoff for P05.CTR.1 : 3 reads## Cutoff for P28.GFP.1 : 9 reads## Cutoff for P28.GFP.2 : 17 reads## Cutoff for P28.MIG.1 : 10 reads## Cutoff for P28.MIG.2 : 5 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 4 reads## Cutoff for P05.GFP.2 : 8 reads## Cutoff for P05.MIG.1 : 3 reads## Cutoff for P05.MIG.2 : 6 reads## Cutoff for P05.CTR.1 : 3 reads## Cutoff for P28.GFP.1 : 9 reads## Cutoff for P28.GFP.2 : 18 reads## Cutoff for P28.MIG.1 : 10 reads## Cutoff for P28.MIG.2 : 5 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 4 reads## Cutoff for P05.GFP.2 : 9 reads## Cutoff for P05.MIG.1 : 3 reads## Cutoff for P05.MIG.2 : 6 reads## Cutoff for P05.CTR.1 : 3 reads## Cutoff for P28.GFP.1 : 9 reads## Cutoff for P28.GFP.2 : 18 reads## Cutoff for P28.MIG.1 : 10 reads## Cutoff for P28.MIG.2 : 5 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 4 reads## Cutoff for P05.GFP.2 : 9 reads## Cutoff for P05.MIG.1 : 3 reads## Cutoff for P05.MIG.2 : 6 reads## Cutoff for P05.CTR.1 : 3 reads## Cutoff for P28.GFP.1 : 9 reads## Cutoff for P28.GFP.2 : 18 reads## Cutoff for P28.MIG.1 : 10 reads## Cutoff for P28.MIG.2 : 5 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 4 reads## Cutoff for P05.GFP.2 : 9 reads## Cutoff for P05.MIG.1 : 3 reads## Cutoff for P05.MIG.2 : 6 reads## Cutoff for P05.CTR.1 : 3 reads## Cutoff for P28.GFP.1 : 9 reads## Cutoff for P28.GFP.2 : 18 reads## Cutoff for P28.MIG.1 : 10 reads## Cutoff for P28.MIG.2 : 6 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 4 reads## Cutoff for P05.GFP.2 : 9 reads## Cutoff for P05.MIG.1 : 3 reads## Cutoff for P05.MIG.2 : 6 reads## Cutoff for P05.CTR.1 : 3 reads## Cutoff for P28.GFP.1 : 10 reads## Cutoff for P28.GFP.2 : 18 reads## Cutoff for P28.MIG.1 : 11 reads## Cutoff for P28.MIG.2 : 6 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 4 reads## Cutoff for P05.GFP.2 : 9 reads## Cutoff for P05.MIG.1 : 3 reads## Cutoff for P05.MIG.2 : 7 reads## Cutoff for P05.CTR.1 : 4 reads## Cutoff for P28.GFP.1 : 10 reads## Cutoff for P28.GFP.2 : 19 reads## Cutoff for P28.MIG.1 : 11 reads## Cutoff for P28.MIG.2 : 6 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 4 reads## Cutoff for P05.GFP.2 : 10 reads## Cutoff for P05.MIG.1 : 3 reads## Cutoff for P05.MIG.2 : 7 reads## Cutoff for P05.CTR.1 : 4 reads## Cutoff for P28.GFP.1 : 10 reads## Cutoff for P28.GFP.2 : 19 reads## Cutoff for P28.MIG.1 : 11 reads## Cutoff for P28.MIG.2 : 6 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 4 reads## Cutoff for P05.GFP.2 : 10 reads## Cutoff for P05.MIG.1 : 4 reads## Cutoff for P05.MIG.2 : 7 reads## Cutoff for P05.CTR.1 : 4 reads## Cutoff for P28.GFP.1 : 10 reads## Cutoff for P28.GFP.2 : 19 reads## Cutoff for P28.MIG.1 : 12 reads## Cutoff for P28.MIG.2 : 6 reads## Cutoff for P28.CTR.1 : 1 reads## Cutoff for P05.GFP.1 : 4 reads## Cutoff for P05.GFP.2 : 10 reads## Cutoff for P05.MIG.1 : 4 reads## Cutoff for P05.MIG.2 : 7 reads## Cutoff for P05.CTR.1 : 4 reads## Cutoff for P28.GFP.1 : 10 reads## Cutoff for P28.GFP.2 : 19 reads## Cutoff for P28.MIG.1 : 12 reads## Cutoff for P28.MIG.2 : 6 reads## Cutoff for P28.CTR.1 : 2 reads## Cutoff for P05.GFP.1 : 5 reads## Cutoff for P05.GFP.2 : 10 reads## Cutoff for P05.MIG.1 : 4 reads## Cutoff for P05.MIG.2 : 7 reads## Cutoff for P05.CTR.1 : 4 reads## Cutoff for P28.GFP.1 : 11 reads## Cutoff for P28.GFP.2 : 20 reads## Cutoff for P28.MIG.1 : 12 reads## Cutoff for P28.MIG.2 : 7 reads## Cutoff for P28.CTR.1 : 2 reads## Cutoff for P05.GFP.1 : 5 reads## Cutoff for P05.GFP.2 : 11 reads## Cutoff for P05.MIG.1 : 4 reads## Cutoff for P05.MIG.2 : 8 reads## Cutoff for P05.CTR.1 : 4 reads## Cutoff for P28.GFP.1 : 11 reads## Cutoff for P28.GFP.2 : 20 reads## Cutoff for P28.MIG.1 : 12 reads## Cutoff for P28.MIG.2 : 7 reads## Cutoff for P28.CTR.1 : 2 reads## Cutoff for P05.GFP.1 : 5 reads## Cutoff for P05.GFP.2 : 11 reads## Cutoff for P05.MIG.1 : 4 reads## Cutoff for P05.MIG.2 : 8 reads## Cutoff for P05.CTR.1 : 4 reads## Cutoff for P28.GFP.1 : 11 reads## Cutoff for P28.GFP.2 : 20 reads## Cutoff for P28.MIG.1 : 13 reads## Cutoff for P28.MIG.2 : 7 reads## Cutoff for P28.CTR.1 : 2 reads## Cutoff for P05.GFP.1 : 5 reads## Cutoff for P05.GFP.2 : 11 reads## Cutoff for P05.MIG.1 : 4 reads## Cutoff for P05.MIG.2 : 8 reads## Cutoff for P05.CTR.1 : 4 reads## Cutoff for P28.GFP.1 : 11 reads## Cutoff for P28.GFP.2 : 21 reads## Cutoff for P28.MIG.1 : 13 reads## Cutoff for P28.MIG.2 : 7 reads## Cutoff for P28.CTR.1 : 2 reads## Cutoff for P05.GFP.1 : 5 reads## Cutoff for P05.GFP.2 : 11 reads## Cutoff for P05.MIG.1 : 4 reads## Cutoff for P05.MIG.2 : 8 reads## Cutoff for P05.CTR.1 : 4 reads## Cutoff for P28.GFP.1 : 12 reads## Cutoff for P28.GFP.2 : 21 reads## Cutoff for P28.MIG.1 : 13 reads## Cutoff for P28.MIG.2 : 7 reads## Cutoff for P28.CTR.1 : 2 reads## Cutoff for P05.GFP.1 : 5 reads## Cutoff for P05.GFP.2 : 11 reads## Cutoff for P05.MIG.1 : 4 reads## Cutoff for P05.MIG.2 : 8 reads## Cutoff for P05.CTR.1 : 5 reads## Cutoff for P28.GFP.1 : 12 reads## Cutoff for P28.GFP.2 : 21 reads## Cutoff for P28.MIG.1 : 14 reads## Cutoff for P28.MIG.2 : 7 reads## Cutoff for P28.CTR.1 : 2 reads## Cutoff for P05.GFP.1 : 5 reads## Cutoff for P05.GFP.2 : 12 reads## Cutoff for P05.MIG.1 : 4 reads## Cutoff for P05.MIG.2 : 9 reads## Cutoff for P05.CTR.1 : 5 reads## Cutoff for P28.GFP.1 : 12 reads## Cutoff for P28.GFP.2 : 21 reads## Cutoff for P28.MIG.1 : 14 reads## Cutoff for P28.MIG.2 : 8 reads## Cutoff for P28.CTR.1 : 2 reads## Cutoff for P05.GFP.1 : 5 reads## Cutoff for P05.GFP.2 : 12 reads## Cutoff for P05.MIG.1 : 5 reads## Cutoff for P05.MIG.2 : 9 reads## Cutoff for P05.CTR.1 : 5 reads## Cutoff for P28.GFP.1 : 13 reads## Cutoff for P28.GFP.2 : 22 reads## Cutoff for P28.MIG.1 : 14 reads## Cutoff for P28.MIG.2 : 8 reads## Cutoff for P28.CTR.1 : 2 reads## Cutoff for P05.GFP.1 : 6 reads## Cutoff for P05.GFP.2 : 12 reads## Cutoff for P05.MIG.1 : 5 reads## Cutoff for P05.MIG.2 : 9 reads## Cutoff for P05.CTR.1 : 5 reads## Cutoff for P28.GFP.1 : 13 reads## Cutoff for P28.GFP.2 : 22 reads## Cutoff for P28.MIG.1 : 15 reads## Cutoff for P28.MIG.2 : 8 reads## Cutoff for P28.CTR.1 : 2 reads## Cutoff for P05.GFP.1 : 6 reads## Cutoff for P05.GFP.2 : 13 reads## Cutoff for P05.MIG.1 : 5 reads## Cutoff for P05.MIG.2 : 9 reads## Cutoff for P05.CTR.1 : 5 reads## Cutoff for P28.GFP.1 : 13 reads## Cutoff for P28.GFP.2 : 22 reads## Cutoff for P28.MIG.1 : 15 reads## Cutoff for P28.MIG.2 : 8 reads## Cutoff for P28.CTR.1 : 2 reads## Cutoff for P05.GFP.1 : 6 reads## Cutoff for P05.GFP.2 : 13 reads## Cutoff for P05.MIG.1 : 5 reads## Cutoff for P05.MIG.2 : 10 reads## Cutoff for P05.CTR.1 : 5 reads## Cutoff for P28.GFP.1 : 14 reads## Cutoff for P28.GFP.2 : 23 reads## Cutoff for P28.MIG.1 : 16 reads## Cutoff for P28.MIG.2 : 9 reads## Cutoff for P28.CTR.1 : 3 reads## Cutoff for P05.GFP.1 : 6 reads## Cutoff for P05.GFP.2 : 13 reads## Cutoff for P05.MIG.1 : 5 reads## Cutoff for P05.MIG.2 : 10 reads## Cutoff for P05.CTR.1 : 5 reads## Cutoff for P28.GFP.1 : 14 reads## Cutoff for P28.GFP.2 : 23 reads## Cutoff for P28.MIG.1 : 16 reads## Cutoff for P28.MIG.2 : 9 reads## Cutoff for P28.CTR.1 : 3 reads## Cutoff for P05.GFP.1 : 6 reads## Cutoff for P05.GFP.2 : 14 reads## Cutoff for P05.MIG.1 : 5 reads## Cutoff for P05.MIG.2 : 10 reads## Cutoff for P05.CTR.1 : 6 reads## Cutoff for P28.GFP.1 : 14 reads## Cutoff for P28.GFP.2 : 23 reads## Cutoff for P28.MIG.1 : 17 reads## Cutoff for P28.MIG.2 : 9 reads## Cutoff for P28.CTR.1 : 3 reads## Cutoff for P05.GFP.1 : 6 reads## Cutoff for P05.GFP.2 : 14 reads## Cutoff for P05.MIG.1 : 5 reads## Cutoff for P05.MIG.2 : 10 reads## Cutoff for P05.CTR.1 : 6 reads## Cutoff for P28.GFP.1 : 15 reads## Cutoff for P28.GFP.2 : 24 reads## Cutoff for P28.MIG.1 : 17 reads## Cutoff for P28.MIG.2 : 9 reads## Cutoff for P28.CTR.1 : 3 reads## Cutoff for P05.GFP.1 : 7 reads## Cutoff for P05.GFP.2 : 14 reads## Cutoff for P05.MIG.1 : 6 reads## Cutoff for P05.MIG.2 : 11 reads## Cutoff for P05.CTR.1 : 6 reads## Cutoff for P28.GFP.1 : 15 reads## Cutoff for P28.GFP.2 : 24 reads## Cutoff for P28.MIG.1 : 18 reads## Cutoff for P28.MIG.2 : 10 reads## Cutoff for P28.CTR.1 : 3 reads## Cutoff for P05.GFP.1 : 7 reads## Cutoff for P05.GFP.2 : 15 reads## Cutoff for P05.MIG.1 : 6 reads## Cutoff for P05.MIG.2 : 11 reads## Cutoff for P05.CTR.1 : 6 reads## Cutoff for P28.GFP.1 : 16 reads## Cutoff for P28.GFP.2 : 25 reads## Cutoff for P28.MIG.1 : 18 reads## Cutoff for P28.MIG.2 : 10 reads## Cutoff for P28.CTR.1 : 3 reads## Cutoff for P05.GFP.1 : 7 reads## Cutoff for P05.GFP.2 : 15 reads## Cutoff for P05.MIG.1 : 6 reads## Cutoff for P05.MIG.2 : 11 reads## Cutoff for P05.CTR.1 : 6 reads## Cutoff for P28.GFP.1 : 16 reads## Cutoff for P28.GFP.2 : 25 reads## Cutoff for P28.MIG.1 : 19 reads## Cutoff for P28.MIG.2 : 11 reads## Cutoff for P28.CTR.1 : 3 reads## Cutoff for P05.GFP.1 : 7 reads## Cutoff for P05.GFP.2 : 16 reads## Cutoff for P05.MIG.1 : 6 reads## Cutoff for P05.MIG.2 : 12 reads## Cutoff for P05.CTR.1 : 7 reads## Cutoff for P28.GFP.1 : 16 reads## Cutoff for P28.GFP.2 : 26 reads## Cutoff for P28.MIG.1 : 19 reads## Cutoff for P28.MIG.2 : 11 reads## Cutoff for P28.CTR.1 : 4 reads## Cutoff for P05.GFP.1 : 7 reads## Cutoff for P05.GFP.2 : 16 reads## Cutoff for P05.MIG.1 : 6 reads## Cutoff for P05.MIG.2 : 12 reads## Cutoff for P05.CTR.1 : 7 reads## Cutoff for P28.GFP.1 : 17 reads## Cutoff for P28.GFP.2 : 26 reads## Cutoff for P28.MIG.1 : 20 reads## Cutoff for P28.MIG.2 : 11 reads## Cutoff for P28.CTR.1 : 4 reads## Cutoff for P05.GFP.1 : 8 reads## Cutoff for P05.GFP.2 : 17 reads## Cutoff for P05.MIG.1 : 7 reads## Cutoff for P05.MIG.2 : 13 reads## Cutoff for P05.CTR.1 : 7 reads## Cutoff for P28.GFP.1 : 18 reads## Cutoff for P28.GFP.2 : 27 reads## Cutoff for P28.MIG.1 : 21 reads## Cutoff for P28.MIG.2 : 12 reads## Cutoff for P28.CTR.1 : 4 reads## Cutoff for P05.GFP.1 : 8 reads## Cutoff for P05.GFP.2 : 17 reads## Cutoff for P05.MIG.1 : 7 reads## Cutoff for P05.MIG.2 : 13 reads## Cutoff for P05.CTR.1 : 7 reads## Cutoff for P28.GFP.1 : 18 reads## Cutoff for P28.GFP.2 : 27 reads## Cutoff for P28.MIG.1 : 22 reads## Cutoff for P28.MIG.2 : 12 reads## Cutoff for P28.CTR.1 : 4 reads## Cutoff for P05.GFP.1 : 8 reads## Cutoff for P05.GFP.2 : 18 reads## Cutoff for P05.MIG.1 : 7 reads## Cutoff for P05.MIG.2 : 14 reads## Cutoff for P05.CTR.1 : 8 reads## Cutoff for P28.GFP.1 : 19 reads## Cutoff for P28.GFP.2 : 28 reads## Cutoff for P28.MIG.1 : 23 reads## Cutoff for P28.MIG.2 : 13 reads## Cutoff for P28.CTR.1 : 5 reads## Cutoff for P05.GFP.1 : 9 reads## Cutoff for P05.GFP.2 : 18 reads## Cutoff for P05.MIG.1 : 7 reads## Cutoff for P05.MIG.2 : 14 reads## Cutoff for P05.CTR.1 : 8 reads## Cutoff for P28.GFP.1 : 20 reads## Cutoff for P28.GFP.2 : 29 reads## Cutoff for P28.MIG.1 : 24 reads## Cutoff for P28.MIG.2 : 13 reads## Cutoff for P28.CTR.1 : 5 reads## Cutoff for P05.GFP.1 : 9 reads## Cutoff for P05.GFP.2 : 19 reads## Cutoff for P05.MIG.1 : 8 reads## Cutoff for P05.MIG.2 : 15 reads## Cutoff for P05.CTR.1 : 8 reads## Cutoff for P28.GFP.1 : 21 reads## Cutoff for P28.GFP.2 : 30 reads## Cutoff for P28.MIG.1 : 25 reads## Cutoff for P28.MIG.2 : 14 reads## Cutoff for P28.CTR.1 : 5 reads## Cutoff for P05.GFP.1 : 9 reads## Cutoff for P05.GFP.2 : 20 reads## Cutoff for P05.MIG.1 : 8 reads## Cutoff for P05.MIG.2 : 16 reads## Cutoff for P05.CTR.1 : 9 reads## Cutoff for P28.GFP.1 : 22 reads## Cutoff for P28.GFP.2 : 30 reads## Cutoff for P28.MIG.1 : 26 reads## Cutoff for P28.MIG.2 : 15 reads## Cutoff for P28.CTR.1 : 6 reads## Cutoff for P05.GFP.1 : 10 reads## Cutoff for P05.GFP.2 : 21 reads## Cutoff for P05.MIG.1 : 9 reads## Cutoff for P05.MIG.2 : 17 reads## Cutoff for P05.CTR.1 : 9 reads## Cutoff for P28.GFP.1 : 23 reads## Cutoff for P28.GFP.2 : 31 reads## Cutoff for P28.MIG.1 : 27 reads## Cutoff for P28.MIG.2 : 16 reads## Cutoff for P28.CTR.1 : 6 reads## Cutoff for P05.GFP.1 : 10 reads## Cutoff for P05.GFP.2 : 22 reads## Cutoff for P05.MIG.1 : 9 reads## Cutoff for P05.MIG.2 : 18 reads## Cutoff for P05.CTR.1 : 10 reads## Cutoff for P28.GFP.1 : 24 reads## Cutoff for P28.GFP.2 : 33 reads## Cutoff for P28.MIG.1 : 29 reads## Cutoff for P28.MIG.2 : 17 reads## Cutoff for P28.CTR.1 : 7 reads## Cutoff for P05.GFP.1 : 11 reads## Cutoff for P05.GFP.2 : 23 reads## Cutoff for P05.MIG.1 : 10 reads## Cutoff for P05.MIG.2 : 19 reads## Cutoff for P05.CTR.1 : 11 reads## Cutoff for P28.GFP.1 : 26 reads## Cutoff for P28.GFP.2 : 34 reads## Cutoff for P28.MIG.1 : 32 reads## Cutoff for P28.MIG.2 : 18 reads## Cutoff for P28.CTR.1 : 7 reads## Cutoff for P05.GFP.1 : 12 reads## Cutoff for P05.GFP.2 : 25 reads## Cutoff for P05.MIG.1 : 11 reads## Cutoff for P05.MIG.2 : 20 reads## Cutoff for P05.CTR.1 : 12 reads## Cutoff for P28.GFP.1 : 29 reads## Cutoff for P28.GFP.2 : 37 reads## Cutoff for P28.MIG.1 : 35 reads## Cutoff for P28.MIG.2 : 20 reads## Cutoff for P28.CTR.1 : 8 reads## Cutoff for P05.GFP.1 : 13 reads## Cutoff for P05.GFP.2 : 27 reads## Cutoff for P05.MIG.1 : 12 reads## Cutoff for P05.MIG.2 : 22 reads## Cutoff for P05.CTR.1 : 13 reads## Cutoff for P28.GFP.1 : 33 reads## Cutoff for P28.GFP.2 : 40 reads## Cutoff for P28.MIG.1 : 41 reads## Cutoff for P28.MIG.2 : 24 reads## Cutoff for P28.CTR.1 : 10 reads## Cutoff for P05.GFP.1 : 15 reads## Cutoff for P05.GFP.2 : 31 reads## Cutoff for P05.MIG.1 : 14 reads## Cutoff for P05.MIG.2 : 26 reads## Cutoff for P05.CTR.1 : 15 readsrownames(table.demux) <- table.demux$sample
table.demux <- table.demux[,2:ncol(table.demux)]
table.demux <- data.frame(t(table.demux))
rownames(table.demux) <- NULL
table.demux## quantile Doublet P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1
## 1 0.50 46325 127 175 144 134 140 187
## 2 0.51 46090 156 205 152 167 158 219
## 3 0.52 45982 167 205 158 181 170 233
## 4 0.53 45982 167 205 158 181 170 233
## 5 0.54 45774 186 234 177 208 190 249
## 6 0.55 45113 216 314 244 285 272 337
## 7 0.56 45002 227 329 242 295 281 351
## 8 0.57 44553 276 356 285 312 323 415
## 9 0.58 44408 293 376 297 323 335 437
## 10 0.59 43900 333 441 354 377 391 433
## 11 0.60 43562 337 493 373 421 426 468
## 12 0.61 43266 374 497 397 450 471 495
## 13 0.62 43266 374 497 397 450 471 495
## 14 0.63 43266 374 497 397 450 471 495
## 15 0.64 43042 389 535 418 447 500 519
## 16 0.65 42430 435 610 459 507 568 596
## 17 0.66 42085 485 621 496 535 611 636
## 18 0.67 41454 543 702 552 597 688 721
## 19 0.68 41364 548 718 553 603 700 730
## 20 0.69 39706 689 912 739 809 766 857
## 21 0.70 38710 798 996 842 880 881 959
## 22 0.71 38710 798 996 842 880 881 959
## 23 0.72 38611 807 1009 840 888 893 971
## 24 0.73 38280 868 1008 879 915 927 1012
## 25 0.74 37899 886 1060 928 962 958 1055
## 26 0.75 37727 900 1082 940 977 976 1076
## 27 0.76 37023 965 1168 1021 1012 1053 1166
## 28 0.77 35964 1048 1260 1129 1124 1168 1220
## 29 0.78 35864 1055 1271 1133 1137 1175 1232
## 30 0.79 35806 1060 1277 1136 1144 1182 1237
## 31 0.80 34516 1180 1399 1262 1270 1249 1381
## 32 0.81 34318 1196 1423 1283 1295 1264 1402
## 33 0.82 34268 1199 1432 1282 1298 1272 1408
## 34 0.83 32832 1346 1581 1429 1445 1428 1516
## 35 0.84 32647 1357 1610 1441 1453 1447 1533
## 36 0.85 32198 1390 1629 1480 1501 1497 1585
## 37 0.86 31882 1406 1685 1508 1522 1540 1621
## 38 0.87 31215 1468 1732 1569 1596 1562 1705
## 39 0.88 30315 1531 1866 1655 1700 1666 1742
## 40 0.89 29804 1579 1898 1710 1735 1733 1791
## 41 0.90 29622 1589 1934 1724 1752 1744 1807
## 42 0.91 28647 1662 2030 1822 1855 1817 1854
## 43 0.92 27908 1728 2079 1908 1952 1892 1930
## 44 0.93 27372 1772 2111 1951 2001 1942 1990
## 45 0.94 26689 1833 2226 2019 2050 1990 2025
## 46 0.95 26146 1884 2250 2062 2089 2041 2094
## 47 0.96 24787 2011 2372 2207 2205 2141 2207
## 48 0.97 23929 2077 2460 2296 2281 2209 2282
## 49 0.98 22459 2215 2543 2427 2397 2315 2419
## 50 0.99 20831 2348 2655 2514 2512 2384 2570
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1 Negative
## 1 126 179 159 155 35
## 2 140 215 167 178 39
## 3 149 229 181 192 39
## 4 149 229 181 192 39
## 5 177 256 203 187 45
## 6 231 291 279 249 55
## 7 241 306 292 263 57
## 8 290 352 346 314 64
## 9 307 372 343 328 67
## 10 374 424 393 389 77
## 11 408 472 424 419 83
## 12 428 508 451 460 89
## 13 428 508 451 460 89
## 14 428 508 451 460 89
## 15 443 544 470 486 93
## 16 506 634 516 522 103
## 17 537 671 546 557 106
## 18 603 662 612 629 123
## 19 617 669 621 640 123
## 20 806 851 811 796 144
## 21 908 968 888 889 167
## 22 908 968 888 889 167
## 23 914 978 905 900 170
## 24 942 1017 932 931 175
## 25 985 1060 963 943 187
## 26 999 1082 974 954 199
## 27 1075 1092 1070 1020 221
## 28 1200 1233 1153 1146 241
## 29 1206 1253 1160 1156 244
## 30 1214 1265 1153 1160 252
## 31 1351 1427 1263 1302 286
## 32 1368 1456 1287 1294 300
## 33 1371 1465 1291 1298 302
## 34 1535 1551 1429 1456 338
## 35 1546 1575 1448 1481 348
## 36 1594 1626 1498 1531 357
## 37 1613 1655 1525 1550 379
## 38 1677 1728 1603 1624 407
## 39 1760 1785 1682 1720 464
## 40 1816 1818 1738 1779 485
## 41 1839 1841 1754 1777 503
## 42 1942 1949 1859 1879 570
## 43 2008 1990 1929 1959 603
## 44 2070 2050 1974 2008 645
## 45 2128 2102 2024 2082 718
## 46 2187 2163 2082 2120 768
## 47 2325 2300 2224 2226 881
## 48 2398 2378 2303 2309 964
## 49 2543 2526 2486 2446 1110
## 50 2683 2691 2626 2584 1488#color.FB <- scales::hue_pal()(10)
par(mfrow=c(3,5))
plot(table.demux$Doublet, type="p", col="#FF6C91", pch=16, ylab="Doublet")
plot(table.demux$Negative, type="p", col="lightgrey", pch=16, ylab="Negative")
plot(rowSums(table.demux[,grep("P",colnames(table.demux))]), type="p", col="#306000", pch=16, ylab="Singlet")
plot(table.demux$quantile, pch=16, ylab="mod-quantile")
plot.new()
plot(table.demux[,2+1], type="p", col=color.FB[1], pch=16, ylab=colnames(table.demux)[2+1])
plot(table.demux[,2+2], type="p", col=color.FB[2], pch=16, ylab=colnames(table.demux)[2+2])
plot(table.demux[,2+3], type="p", col=color.FB[3], pch=16, ylab=colnames(table.demux)[2+3])
plot(table.demux[,2+4], type="p", col=color.FB[4], pch=16, ylab=colnames(table.demux)[2+4])
plot(table.demux[,2+5], type="p", col=color.FB[5], pch=16, ylab=colnames(table.demux)[2+5])
plot(table.demux[,2+6], type="p", col=color.FB[6], pch=16, ylab=colnames(table.demux)[2+6])
plot(table.demux[,2+7], type="p", col=color.FB[7], pch=16, ylab=colnames(table.demux)[2+7])
plot(table.demux[,2+8], type="p", col=color.FB[8], pch=16, ylab=colnames(table.demux)[2+8])
plot(table.demux[,2+9], type="p", col=color.FB[9], pch=16, ylab=colnames(table.demux)[2+9])
plot(table.demux[,2+10], type="p", col=color.FB[10], pch=16, ylab=colnames(table.demux)[2+10])
q0.9900 ~ 0.9999
#color.FB <- scales::hue_pal()(10)
par(mfrow=c(3,5))
plot(table.demux.s$Doublet, type="p", col="#FF6C91", pch=16, ylab="Doublet")
plot(table.demux.s$Negative, type="p", col="lightgrey", pch=16, ylab="Negative")
plot(rowSums(table.demux.s[,grep("P",colnames(table.demux.s))]), type="p", col="#306000", pch=16, ylab="Singlet")
plot(table.demux.s$quantile, pch=16, ylab="mod-quantile")
plot.new()
plot(table.demux.s[,2+1], type="p", col=color.FB[1], pch=16, ylab=colnames(table.demux.s)[2+1])
plot(table.demux.s[,2+2], type="p", col=color.FB[2], pch=16, ylab=colnames(table.demux.s)[2+2])
plot(table.demux.s[,2+3], type="p", col=color.FB[3], pch=16, ylab=colnames(table.demux.s)[2+3])
plot(table.demux.s[,2+4], type="p", col=color.FB[4], pch=16, ylab=colnames(table.demux.s)[2+4])
plot(table.demux.s[,2+5], type="p", col=color.FB[5], pch=16, ylab=colnames(table.demux.s)[2+5])
plot(table.demux.s[,2+6], type="p", col=color.FB[6], pch=16, ylab=colnames(table.demux.s)[2+6])
plot(table.demux.s[,2+7], type="p", col=color.FB[7], pch=16, ylab=colnames(table.demux.s)[2+7])
plot(table.demux.s[,2+8], type="p", col=color.FB[8], pch=16, ylab=colnames(table.demux.s)[2+8])
plot(table.demux.s[,2+9], type="p", col=color.FB[9], pch=16, ylab=colnames(table.demux.s)[2+9])
plot(table.demux.s[,2+10], type="p", col=color.FB[10], pch=16, ylab=colnames(table.demux.s)[2+10])
demultiplexing
(tags1/2/6-10 q0.998; tags3-5 q0.99)
# Demultiplex cells based on HTO enrichment
FB.seur <- HTODemux(FB.seur, assay = 'HTO', positive.quantile = 0.994)## Cutoff for P28.GFP.1 : 36 reads## Cutoff for P28.GFP.2 : 43 reads## Cutoff for P28.MIG.1 : 45 reads## Cutoff for P28.MIG.2 : 26 reads## Cutoff for P28.CTR.1 : 12 reads## Cutoff for P05.GFP.1 : 17 reads## Cutoff for P05.GFP.2 : 34 reads## Cutoff for P05.MIG.1 : 15 reads## Cutoff for P05.MIG.2 : 29 reads## Cutoff for P05.CTR.1 : 17 readstable(FB.seur$HTO_classification.global)##
## Doublet Negative Singlet
## 19627 1851 26408table(FB.seur$hash.ID)[c("Doublet","Negative",rownames(FB))]##
## Doublet Negative P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1
## 19627 1851 2461 2741 2559 2613 2374 2672
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1
## 2774 2782 2745 2687# Demultiplex cells based on HTO enrichment
FB.seur <- HTODemux(FB.seur, assay = 'HTO', positive.quantile = 0.996)## Cutoff for P28.GFP.1 : 38 reads## Cutoff for P28.GFP.2 : 45 reads## Cutoff for P28.MIG.1 : 48 reads## Cutoff for P28.MIG.2 : 28 reads## Cutoff for P28.CTR.1 : 13 reads## Cutoff for P05.GFP.1 : 18 reads## Cutoff for P05.GFP.2 : 37 reads## Cutoff for P05.MIG.1 : 16 reads## Cutoff for P05.MIG.2 : 31 reads## Cutoff for P05.CTR.1 : 18 readstable(FB.seur$HTO_classification.global)##
## Doublet Negative Singlet
## 18932 2130 26824table(FB.seur$hash.ID)[c("Doublet","Negative",rownames(FB))]##
## Doublet Negative P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1
## 18932 2130 2524 2779 2564 2645 2363 2739
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1
## 2807 2851 2810 2742# Demultiplex cells based on HTO enrichment
FB.seur <- HTODemux(FB.seur, assay = 'HTO', positive.quantile = 0.998)## Cutoff for P28.GFP.1 : 43 reads## Cutoff for P28.GFP.2 : 48 reads## Cutoff for P28.MIG.1 : 54 reads## Cutoff for P28.MIG.2 : 32 reads## Cutoff for P28.CTR.1 : 15 reads## Cutoff for P05.GFP.1 : 20 reads## Cutoff for P05.GFP.2 : 41 reads## Cutoff for P05.MIG.1 : 18 reads## Cutoff for P05.MIG.2 : 35 reads## Cutoff for P05.CTR.1 : 20 readstable(FB.seur$HTO_classification.global)##
## Doublet Negative Singlet
## 17761 2819 27306table(FB.seur$hash.ID)[c("Doublet","Negative",rownames(FB))]##
## Doublet Negative P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1
## 17761 2819 2629 2865 2495 2643 2284 2839
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1
## 2849 2943 2909 2850# Demultiplex cells based on HTO enrichment
FB.seur <- HTODemux(FB.seur, assay = 'HTO', positive.quantile = 0.99)## Cutoff for P28.GFP.1 : 33 reads## Cutoff for P28.GFP.2 : 40 reads## Cutoff for P28.MIG.1 : 41 reads## Cutoff for P28.MIG.2 : 24 reads## Cutoff for P28.CTR.1 : 10 reads## Cutoff for P05.GFP.1 : 15 reads## Cutoff for P05.GFP.2 : 31 reads## Cutoff for P05.MIG.1 : 14 reads## Cutoff for P05.MIG.2 : 26 reads## Cutoff for P05.CTR.1 : 15 readstable(FB.seur$HTO_classification.global)##
## Doublet Negative Singlet
## 20831 1488 25567table(FB.seur$hash.ID)[c("Doublet","Negative",rownames(FB))]##
## Doublet Negative P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1
## 20831 1488 2348 2655 2514 2512 2384 2570
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1
## 2683 2691 2626 2584# tags q0.99
#cutoff.FB <- c(33,40,41,24,10,
# 15,31,14,26,15)
# tags1/2/6-10 q0.996; tags3-5 q0.99
#cutoff.FB <- c(38,45,41,24,10,
# 18,37,16,31,18)
# tags1/2/6-10 q0.998; tags3-5 q0.99
cutoff.FB <- c(43,48,41,24,10,
20,41,18,35,20)
names(cutoff.FB) <- rownames(FB.seur)
cutoff.FB## P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1 P05.GFP.2 P05.MIG.1
## 43 48 41 24 10 20 41 18
## P05.MIG.2 P05.CTR.1
## 35 20par(mfrow=c(2,5))
#plot(sort(rowSums(t(FB))),
# xlab="reordered cells",
# ylab="UMI counts", log="xy",
# main="sum")
plot(sort(t(FB)[,1]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[1])
abline(h=cutoff.FB[1],col = color.FB[1])
plot(sort(t(FB)[,2]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[2])
abline(h=cutoff.FB[2],col = color.FB[2])
plot(sort(t(FB)[,3]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[3])
abline(h=cutoff.FB[3],col = color.FB[3])
plot(sort(t(FB)[,4]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[4])
abline(h=cutoff.FB[4],col = color.FB[4])
plot(sort(t(FB)[,5]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[5])
abline(h=cutoff.FB[5],col = color.FB[5])
plot(sort(t(FB)[,6]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[6])
abline(h=cutoff.FB[6],col = color.FB[6])
plot(sort(t(FB)[,7]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[7])
abline(h=cutoff.FB[7],col = color.FB[7])
plot(sort(t(FB)[,8]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[8])
abline(h=cutoff.FB[8],col = color.FB[8])
plot(sort(t(FB)[,9]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[9])
abline(h=cutoff.FB[9],col = color.FB[9])
plot(sort(t(FB)[,10]),
xlab="reordered cells",
ylab="UMI counts", log="xy",
main=rownames(FB)[10])
abline(h=cutoff.FB[10],col = color.FB[10])
custom.Demux <- function(FBobj,FBcutoff){
# define decision matrix
dd <- FBobj@assays[['HTO']]@counts
dd <- apply(dd, 2, function(x){
x > FBcutoff
})
x1 <- colSums(dd)
x1[x1>1] <- "Doublet"
x1[x1==0] <- "Negative"
x1[x1==1] <- "Singlet"
x2 <- x1
bc.slt <- names(x1)[x1=="Singlet"]
for(bc in bc.slt){
x2[bc] <- rownames(dd)[dd[,bc]]
}
FBdf <- data.frame(HTO_classification.global=factor(x1,
levels = c("Doublet","Negative","Singlet")),
hash.ID=factor(x2,
levels = c("Doublet","Negative",rownames(dd))))
return(FBdf)
}
df.FB <- custom.Demux(FB.seur,cutoff.FB)dim(df.FB)## [1] 47886 2df.FB[1:10,]## HTO_classification.global hash.ID
## AAACCCAAGAAATCCA-1 Singlet P05.GFP.1
## AAACCCAAGAATCGTA-1 Singlet P28.GFP.1
## AAACCCAAGAGCCTGA-1 Doublet Doublet
## AAACCCAAGATGCTAA-1 Doublet Doublet
## AAACCCAAGCACCTGC-1 Doublet Doublet
## AAACCCAAGCGATCGA-1 Singlet P05.MIG.2
## AAACCCAAGGTAGCCA-1 Singlet P28.GFP.1
## AAACCCAAGGTTGAGC-1 Singlet P28.GFP.1
## AAACCCAAGTAAACGT-1 Doublet Doublet
## AAACCCAAGTAAGCAT-1 Doublet Doublettable(df.FB)## hash.ID
## HTO_classification.global Doublet Negative P28.GFP.1 P28.GFP.2 P28.MIG.1
## Doublet 18820 0 0 0 0
## Negative 0 1966 0 0 0
## Singlet 0 0 2496 2793 2698
## hash.ID
## HTO_classification.global P28.MIG.2 P28.CTR.1 P05.GFP.1 P05.GFP.2 P05.MIG.1
## Doublet 0 0 0 0 0
## Negative 0 0 0 0 0
## Singlet 2732 2594 2721 2725 2836
## hash.ID
## HTO_classification.global P05.MIG.2 P05.CTR.1
## Doublet 0 0
## Negative 0 0
## Singlet 2785 2720FB.seur@meta.data[,c("HTO_classification.global","hash.ID")] <- df.FB
FB.seur@meta.data[1:4,]## orig.ident nCount_RNA nFeature_RNA nCount_HTO nFeature_HTO
## AAACCCAAGAAATCCA-1 plus_SIMPLE 221 10 221 10
## AAACCCAAGAATCGTA-1 plus_SIMPLE 80 8 80 8
## AAACCCAAGAGCCTGA-1 plus_SIMPLE 670 10 670 10
## AAACCCAAGATGCTAA-1 plus_SIMPLE 142 9 142 9
## HTO_maxID HTO_secondID HTO_margin HTO_classification
## AAACCCAAGAAATCCA-1 P05.GFP.1 P05.GFP.2 2.0828322 P05.GFP.1
## AAACCCAAGAATCGTA-1 P28.GFP.1 P28.MIG.1 1.0297322 P28.GFP.1
## AAACCCAAGAGCCTGA-1 P28.GFP.1 P05.MIG.2 0.1030854 P05.MIG.2_P28.GFP.1
## AAACCCAAGATGCTAA-1 P05.CTR.1 P28.MIG.2 0.1557043 P05.CTR.1_P28.MIG.2
## HTO_classification.global hash.ID
## AAACCCAAGAAATCCA-1 Singlet P05.GFP.1
## AAACCCAAGAATCGTA-1 Singlet P28.GFP.1
## AAACCCAAGAGCCTGA-1 Doublet Doublet
## AAACCCAAGATGCTAA-1 Doublet Doublet## if all using the same q, run this one
#FB.seur$HTO_classification.global <- factor(as.character(FB.seur$HTO_classification.global),
# levels = c("Doublet","Negative","Singlet"))
#FB.seur$hash.ID <- factor(as.character(FB.seur$hash.ID),
# levels = c("Doublet","Negative",rownames(FB.seur)))sl_stat <- table(FB.seur$HTO_classification.global)
barplot(sl_stat,ylim = c(0,32500),col = c("#FF6C91","lightgrey","#7CAE00"),main = "Feature Barcode statistics")
text(x=1:3*1.2-0.5,y=sl_stat+1680,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)
RidgePlot
with ridge plots, raw
FB.seur@active.ident <- factor(FB.seur$HTO_maxID,levels=rownames(FB.seur))
RidgePlot(FB.seur, assay = 'HTO', features = rownames(FB.seur[['HTO']]),same.y.lims= TRUE, ncol=5,
cols = color.FB) 
with ridge plots, filtered
#FB.seur@meta.data$hash.ID.sort <- factor(as.character(FB.seur@meta.data$hash.ID),levels = c("Doublet","Negative",levels(FB.seur@active.ident)))
RidgePlot(FB.seur, assay = 'HTO', features = rownames(FB.seur[['HTO']]), ncol=5,same.y.lims= TRUE,
cols = c("#FF6C91","lightgrey",color.FB),group.by = "hash.ID") 
tSNE for HTOs
# first, remove negative cells from th object
#FB.seur.subset <- FB.seur
# Calculate a tSNE embedding of the HTO data
#DefaultAssay(FB.seur.subset) <- "HTO"
#FB.seur.subset <- ScaleData(FB.seur.subset, features = rownames(FB.seur.subset), verbose = FALSE)
#FB.seur.subset <- RunPCA(FB.seur.subset, features = rownames(FB.seur.subset), approx = FALSE)
#FB.seur.subset <- RunTSNE(FB.seur.subset, dims = 1:10, perplexity = 500, check_duplicates = FALSE)
#saveRDS(FB.seur.subset, "FB0811.seur.subset.rds")
FB.seur.subset <- readRDS("FB0811.seur.subset.rds")multiplot(
DimPlot(FB.seur.subset, order = rev(rownames(FB.seur)),
cols = color.FB) + labs(title = "FB Max_ID"),
DimPlot(FB.seur.subset, order = rev(c("Negative",rownames(FB.seur),"Doublet")),
group.by = 'hash.ID', label = F,
cols = c("lightgrey",color.FB,"#FF6C91")) +
labs(title = "FB Filtered_ID") + theme(plot.title = element_text(hjust = 0)) ,
cols = 1)
stat
FB.seur$HTO_classification.global <- factor(as.character(FB.seur$HTO_classification.global),
levels = c("Doublet","Negative","Singlet"))
Idents(FB.seur) <- "HTO_classification.global"
VlnPlot(FB.seur, features = "nCount_RNA", pt.size = 0, log = TRUE, cols = c("#F8766D","lightgrey","#00BA38"))
VlnPlot(FB.seur, features = "nCount_RNA", pt.size = 0, log = TRUE, cols = c("#F8766D","lightgrey","#00BA38")) + geom_jitter(alpha=0.015)
table(FB.seur@meta.data$hash.ID)##
## Doublet Negative P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1
## 18820 1966 2496 2793 2698 2732 2594 2721
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1
## 2725 2836 2785 2720par(mar=c(6,4,2,3))
sl_stat <- table(FB.seur$hash.ID)
barplot(sl_stat,ylim = c(0,23800),
col = c("#FF6C91","lightgrey",color.FB),
main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:12*1.2-0.48 ,levels(FB.seur@meta.data$hash.ID.sort), las=3, cex.axis=0.85)
text(x=1:12*1.2-0.45,y=sl_stat+1575,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)
GEX
mainly follow https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
GEX.seur <- CreateSeuratObject(counts = GEX,
min.cells = 3,
min.features = 200,
project = "plus_SIMPLE")
GEX.seur## An object of class Seurat
## 18371 features across 47886 samples within 1 assay
## Active assay: RNA (18371 features, 0 variable features)filtering
MT genes
grep("^mt-",rownames(GEX),value = T)## [1] "mt-Nd1" "mt-Nd2" "mt-Co1" "mt-Co2" "mt-Atp8" "mt-Atp6" "mt-Co3"
## [8] "mt-Nd3" "mt-Nd4l" "mt-Nd4" "mt-Nd5" "mt-Nd6" "mt-Cytb"MT_gene <- grep("^mt-",rownames(GEX.seur),value = T)
GEX.seur[["percent.mt"]] <- PercentageFeatureSet(GEX.seur, features = MT_gene)
RB_gene <- grep("^Rps|^Rpl",rownames(GEX.seur),value = T)
GEX.seur[["percent.rb"]] <- PercentageFeatureSet(GEX.seur, features = RB_gene)
# Visualize QC metrics as a violin plot
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0.01)
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0)
plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt")
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb")
plota + plotb + plotc
add FB info
GEX.seur$FB.info <- FB.seur@meta.data[rownames(GEX.seur@meta.data),"hash.ID"]
#head(GEX.seur@meta.data)table(GEX.seur$FB.info)##
## Doublet Negative P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1
## 18820 1966 2496 2793 2698 2732 2594 2721
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1
## 2725 2836 2785 2720VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = c("#FF6C91","lightgrey",color.FB),
pt.size = 0.1, split.by = "FB.info")## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##
## This message will be shown once per session.
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = c("#FF6C91","lightgrey",color.FB),
pt.size = 0, split.by = "FB.info")
GEX.seur <- subset(GEX.seur, subset = FB.info!="Doublet" & FB.info!="Negative")
GEX.seur## An object of class Seurat
## 18371 features across 27100 samples within 1 assay
## Active assay: RNA (18371 features, 0 variable features)VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0.1)
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0)
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = c("#FF6C91","lightgrey",color.FB),
pt.size = 0.1, split.by = "FB.info")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = c("#FF6C91","lightgrey",color.FB),
pt.size = 0, split.by = "FB.info")
plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt")
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb")
plota + plotb + plotc
filtering
GEX.seur <- subset(GEX.seur, subset = percent.mt < 10 & nFeature_RNA > 800 & nFeature_RNA < 3200 & nCount_RNA < 10000)
GEX.seur## An object of class Seurat
## 18371 features across 25338 samples within 1 assay
## Active assay: RNA (18371 features, 0 variable features)filtered obj
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0.1)
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0)
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = color.FBraw,
pt.size = 0.1, split.by = "FB.info")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
#ncol = 3,
ncol = 2,
cols = color.FBraw,
pt.size = 0, split.by = "FB.info")
plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt")
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb")
plota + plotb + plotc
par(mar=c(6,4,2,3))
sl_stat <- table(FB.seur$hash.ID)
barplot(sl_stat,ylim = c(0,23800),
col = c("#FF6C91","lightgrey",color.FB),
main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:12*1.2-0.48 ,levels(FB.seur@meta.data$hash.ID), las=3, cex.axis=0.85)
text(x=1:12*1.2-0.45,y=sl_stat+1575,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)
par(mar=c(6,4,2,3))
sl_stat <- table(GEX.seur$FB.info)
barplot(sl_stat,ylim = c(0,4800),
col = c("#FF6C91","lightgrey",color.FB),
main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:12*1.2-0.48 ,levels(FB.seur@meta.data$hash.ID), las=3, cex.axis=0.85)
text(x=1:12*1.2-0.45,y=sl_stat+345,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)
check cellcycle
s.genes <- Hmisc::capitalize(tolower(cc.genes.updated.2019$s.genes))
g2m.genes <- Hmisc::capitalize(tolower(cc.genes.updated.2019$g2m.genes))
GEX.seur <- CellCycleScoring(GEX.seur,
s.features = s.genes,
g2m.features = g2m.genes)VlnPlot(GEX.seur, features = c("S.Score", "G2M.Score"), group.by = "FB.info",
ncol = 2, pt.size = 0.1)
GEX.seur.cc <- GEX.seur
GEX.seur.cc <- NormalizeData(GEX.seur.cc)
GEX.seur.cc <- FindVariableFeatures(GEX.seur.cc)
GEX.seur.cc <- ScaleData(GEX.seur.cc)
Idents(GEX.seur.cc) <- "Phase"
RidgePlot(GEX.seur.cc, features = c("Pcna", "Top2a", "Mcm6", "Mki67"), ncol = 2)
GEX.seur.cc <- RunPCA(GEX.seur.cc, features = c(s.genes, g2m.genes))
DimPlot(GEX.seur.cc, reduction = 'pca')
seems there’re partly cycling !
Markers and Clusters
Normalizing
GEX.seur <- NormalizeData(GEX.seur, normalization.method = "LogNormalize", scale.factor = 10000)GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1500)
# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2
head(VariableFeatures(GEX.seur), 200)## [1] "Cd74" "Spp1" "H2-Aa" "H2-Eb1"
## [5] "H2-Ab1" "Ccl4" "Cxcl10" "Ccl5"
## [9] "Pf4" "Ccl3" "Spp1-EGFP" "S100a6"
## [13] "Ccl12" "Hist1h1b" "Crip1" "S100a4"
## [17] "Mki67" "Cxcl9" "Ccl2" "Cd209a"
## [21] "Ifi27l2a" "Cxcl13" "Rsad2" "Top2a"
## [25] "Apoe" "Fn1" "Xist" "Lgals3"
## [29] "Ccr2" "Fabp5" "Plac8" "Ifitm3"
## [33] "Cst7" "Lpl" "Lyve1" "Hist1h2ae"
## [37] "Lyz2" "Cenpf" "S100a11" "Mgl2"
## [41] "F13a1" "Hmgb2" "Hist1h2ap" "Ms4a7"
## [45] "Meg3" "Snhg11" "Ch25h" "Ccl7"
## [49] "Ifit3" "Isg15" "Cd163" "Ifi203"
## [53] "Prc1" "Ifit2" "AA467197" "Vim"
## [57] "Cd209f" "Stmn1" "Iigp1" "Gpnmb"
## [61] "Ear2" "Mrc1" "Syt1" "Pclaf"
## [65] "Wfdc17" "Cd52" "Ifit1" "Napsa"
## [69] "Ube2c" "Egr1" "Lgals1" "Hmmr"
## [73] "Mmp12" "S100a10" "Plbd1" "Ahnak"
## [77] "Ccl8" "Nusap1" "Il1b" "Apoc1"
## [81] "Nkg7" "Cenpa" "Cd72" "Arg1"
## [85] "Cd69" "Cybb" "Csf1" "Hmox1"
## [89] "Gria2" "Ccnd2" "Nrxn3" "Cd3d"
## [93] "Atp1b1" "Aspm" "Fxyd5" "Cenpe"
## [97] "Ccl6" "Ifitm6" "Birc5" "Apoc2"
## [101] "Gm26737" "Rian" "Cd9" "Serpina3g"
## [105] "Epcam" "Ctsd" "Apod" "Ccr7"
## [109] "Mndal" "Cdca8" "Cd3g" "Cadps"
## [113] "Pbld1" "H2-K1" "Ptprcap" "Igf1"
## [117] "H2afz" "Ifi204" "Ms4a4b" "Hist1h3c"
## [121] "Tpx2" "Spn" "Ptgds" "Olfm1"
## [125] "Ndn" "Lgals3bp" "Slamf7" "Trbc2"
## [129] "Gdf15" "Il12b" "Fth1" "H2-D1"
## [133] "Mif" "Bst2" "Klrd1" "Igfbp5"
## [137] "Clec4b1" "Slpi" "Stmn3" "Gad1"
## [141] "Smc2" "Ly6a" "Saa3" "Ctsb"
## [145] "Clec10a" "Nfasc" "Ptprn" "Kif11"
## [149] "Kctd16" "Slfn1" "B230312C02Rik" "Ifitm2"
## [153] "Stmn2" "Acod1" "Ttr" "Gzmb"
## [157] "Kctd14" "Esco2" "Ifi205" "Slfn5"
## [161] "Igfbp2" "Hp" "Cd3e" "H2afx"
## [165] "Rcan1" "Cdkn1a" "Knl1" "Smc4"
## [169] "Chl1" "Ifi211" "Ace" "Il2rb"
## [173] "Gimap3" "Ccnb1" "Rrm2" "Ftl1"
## [177] "Hist1h1e" "Ank3" "Cd24a" "Traf1"
## [181] "Adarb2" "B830012L14Rik" "Pbk" "Emb"
## [185] "Tmsb10" "Gad2" "Hba-a2" "Dab2"
## [189] "Cd209g" "B2m" "Gimap7" "Hba-a1"
## [193] "Ms4a4c" "Trac" "Cldn11" "Cd36"
## [197] "Cd63" "Jaml" "Runx3" "Cd8b1"GEX.seur <- ScaleData(GEX.seur, features = rownames(GEX.seur))## Centering and scaling data matrixPCA
# exclude MT genes and more
#DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^Gm|^Hist",rownames(GEX.seur),value = T)
DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^AC|^AA|^AW|^BC|^Gm|^Hist|Rik$|-ps",
rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
c(MT_gene,
DIG,
CC_gene) )
GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur))## PC_ 1
## Positive: Cst7, Ifitm3, Ifi27l2a, Ifi204, Isg15, B2m, Cd72, Lgals3bp, Fth1, Ch25h
## Cd52, Bst2, H2-K1, Ifi207, Ifit2, Oasl2, Ctsb, Fxyd5, Ifit3, H2-D1
## Zbp1, Rsad2, Ccl3, Irf7, Ctsd, Ctsz, Slfn5, Tspo, Ccl2, Cd9
## Negative: Crybb1, Ccr5, Sox4, Fchsd2, Mrc1, Igfbp4, Rgs7bp, Hmgb1, Cryba4, Cbfa2t3
## Il7r, Gbp7, Gas6, Atrx, Stab1, Ccr1, Zbtb20, Khdrbs3, Tmem52, Upk1b
## Slc12a2, Stox2, Bod1l, Lsp1, H1f0, Clec4a1, Smc3, Lcorl, Ednrb, Cask
## PC_ 2
## Positive: Prc1, Pclaf, Pbk, Knl1, Kif15, Esco2, Aspm, Racgap1, Spc24, Lockd
## Mis18bp1, Stmn1, Smc2, Bub1b, Neil3, Ccna2, Sgo2a, Ccnb1, Ccnf, Plk1
## H2afx, Spc25, Cit, Tk1, Ncapg, Shcbp1, Diaph3, Kif4, Trim59, Kif22
## Negative: Ctsd, Tyrobp, Ctsb, Timp2, Apoe, Ctsl, Cd63, B2m, Abca1, H2-D1
## H2-K1, Cd9, Ctsz, Mpeg1, Fth1, Lgals3bp, Grn, Npc2, Ftl1, Lag3
## Lyz2, Pmp22, Xist, H2-T23, Slfn5, Ly86, Cst7, Ifi204, Ccl9, Abcg1
## PC_ 3
## Positive: Ctsd, Cd9, Ctsl, Ctsb, Ctsz, Syngr1, Cst7, Cd63, Serpine2, Lpl
## Tyrobp, Grn, Timp2, Ly86, Ccl3, Fabp5, Mif, Baiap2l2, Pld3, Atp6v0c
## Hif1a, Csf1, Spp1-EGFP, Pmp22, Gnas, Cadm1, Myo1e, Spp1, Fam20c, Cpd
## Negative: H2-Aa, H2-Eb1, S100a6, Ifi203, Plbd1, H2-Ab1, S100a11, Cd74, S100a4, Mndal
## Crip1, Ccr2, Cybb, Vim, Napsa, Ahnak, Clec12a, Slamf7, Itga4, S100a10
## Emb, Mgl2, Ms4a4c, Clec10a, Ciita, Ifi211, Wfdc17, Ifitm2, Ifit3b, Ifit3
## PC_ 4
## Positive: Lgals3, Crip1, Lgals1, Vim, Fxyd5, H2-Ab1, Ftl1, H2-Eb1, Anxa2, H2-Aa
## S100a6, S100a11, S100a4, Cd74, Igf1, Plbd1, Iqgap1, Ifitm2, Ccr2, Anxa5
## Alcam, Ahnak, Wfdc17, Napsa, Lpl, Itga4, Fabp5, Bhlhe40, Emp3, Clec12a
## Negative: Ifit3, Ifit2, Ifit1, Rsad2, Iigp1, Ifit3b, Isg15, Ifi204, Cxcl10, Slfn5
## Usp18, Parp14, Oasl2, Phf11d, Ifi206, Ifi209, Herc6, Ccl12, Stat1, Ifi213
## Oasl1, Fgl2, Rtp4, Ifi208, Slfn8, Ifi211, Trim30a, Irf7, Gbp2, Tor3a
## PC_ 5
## Positive: H2-Ab1, H2-Aa, H2-Eb1, S100a11, S100a6, Cd74, Klrd1, Plbd1, Ccr2, S100a4
## Cd9, Napsa, Vim, Crip1, Slamf7, S100a10, H2-DMa, Il7r, Mcemp1, Csmd3
## Alcam, Ciita, Ctsd, Cytip, H2-DMb1, Emb, Olfm1, Anxa2, Cd34, Cd209a
## Negative: Pf4, Ms4a7, Mrc1, F13a1, Cbr2, Dab2, Pla2g7, Cd163, Lyve1, Apoe
## Cp, Ednrb, Gas6, Folr2, Clec4n, Colec12, Igfbp4, Cd36, Pltp, Stab1
## Igf1, Cd38, Itsn1, Msr1, Clec4a1, Gpx3, Myo5a, Clec10a, Adam33, Mgl2length(VariableFeatures(GEX.seur))## [1] 1229head(VariableFeatures(GEX.seur),300)## [1] "Cd74" "Spp1" "H2-Aa" "H2-Eb1" "H2-Ab1" "Ccl4"
## [7] "Cxcl10" "Ccl5" "Pf4" "Ccl3" "Spp1-EGFP" "S100a6"
## [13] "Ccl12" "Crip1" "S100a4" "Cxcl9" "Ccl2" "Cd209a"
## [19] "Ifi27l2a" "Cxcl13" "Rsad2" "Apoe" "Fn1" "Xist"
## [25] "Lgals3" "Ccr2" "Fabp5" "Plac8" "Ifitm3" "Cst7"
## [31] "Lpl" "Lyve1" "Lyz2" "S100a11" "Mgl2" "F13a1"
## [37] "Ms4a7" "Meg3" "Snhg11" "Ch25h" "Ccl7" "Ifit3"
## [43] "Isg15" "Cd163" "Ifi203" "Prc1" "Ifit2" "Vim"
## [49] "Cd209f" "Stmn1" "Iigp1" "Gpnmb" "Ear2" "Mrc1"
## [55] "Syt1" "Pclaf" "Wfdc17" "Cd52" "Ifit1" "Napsa"
## [61] "Egr1" "Lgals1" "Mmp12" "S100a10" "Plbd1" "Ahnak"
## [67] "Ccl8" "Il1b" "Apoc1" "Nkg7" "Cd72" "Arg1"
## [73] "Cd69" "Cybb" "Csf1" "Hmox1" "Gria2" "Ccnd2"
## [79] "Nrxn3" "Cd3d" "Atp1b1" "Aspm" "Fxyd5" "Ccl6"
## [85] "Ifitm6" "Apoc2" "Rian" "Cd9" "Serpina3g" "Epcam"
## [91] "Ctsd" "Apod" "Ccr7" "Mndal" "Cd3g" "Cadps"
## [97] "Pbld1" "H2-K1" "Ptprcap" "Igf1" "H2afz" "Ifi204"
## [103] "Ms4a4b" "Spn" "Ptgds" "Olfm1" "Ndn" "Lgals3bp"
## [109] "Slamf7" "Gdf15" "Il12b" "Fth1" "H2-D1" "Mif"
## [115] "Bst2" "Klrd1" "Igfbp5" "Clec4b1" "Slpi" "Stmn3"
## [121] "Gad1" "Smc2" "Ly6a" "Saa3" "Ctsb" "Clec10a"
## [127] "Nfasc" "Ptprn" "Kctd16" "Slfn1" "Ifitm2" "Stmn2"
## [133] "Acod1" "Ttr" "Gzmb" "Kctd14" "Esco2" "Ifi205"
## [139] "Slfn5" "Igfbp2" "Hp" "Cd3e" "H2afx" "Rcan1"
## [145] "Cdkn1a" "Knl1" "Chl1" "Ifi211" "Ace" "Il2rb"
## [151] "Gimap3" "Ccnb1" "Ftl1" "Ank3" "Cd24a" "Adarb2"
## [157] "Pbk" "Emb" "Tmsb10" "Gad2" "Dab2" "Cd209g"
## [163] "B2m" "Gimap7" "Ms4a4c" "Cldn11" "Cd36" "Cd63"
## [169] "Jaml" "Runx3" "Cd8b1" "Kit" "Ciita" "Ncam1"
## [175] "Adgre5" "Ifi207" "Grik1" "Cd38" "Lilrb4a" "Itgal"
## [181] "Cd5l" "Nr4a3" "Clstn3" "Kif15" "Plp1" "Mis18bp1"
## [187] "Ly6h" "Ctsz" "Slc6a1" "Krt80" "Gbp2" "Grm5"
## [193] "Lsp1" "Clec12a" "Upb1" "Enpp2" "Dqx1" "Ifi44"
## [199] "Mt1" "Ccl9" "Syngr1" "Cbr2" "Sh2d2a" "Tnr"
## [205] "H2-Q7" "Itga4" "Snca" "Acp5" "Id2" "Cd226"
## [211] "Drd4" "Mcemp1" "Ifit3b" "Ccr1" "Mcf2l" "Nsg1"
## [217] "Pla2g2d" "Tspo" "Cxcl2" "Adgre4" "Ccl24" "Irf7"
## [223] "Marco" "C4b" "Cspg4" "Bex2" "Cd93" "Bub1b"
## [229] "Map1b" "Racgap1" "Nap1l5" "Sgo2a" "Nap1l3" "Dennd5b"
## [235] "Wdcp" "Ctsl" "Adamdec1" "Gapdh" "Opcml" "Atf3"
## [241] "Tc2n" "Oasl2" "Anxa1" "Ctla2a" "Cmpk2" "Gnas"
## [247] "Klrb1c" "Rbms3" "Pcdh17" "Aldh1a3" "Scn2a" "Vstm2l"
## [253] "Cstb" "Adcy1" "Clec4n" "Sst" "Spc24" "Peg3"
## [259] "Gap43" "Baiap2l2" "Brinp2" "Ppp1r14a" "Rnase2a" "Ctsw"
## [265] "Sct" "Ifitm1" "Fpr1" "Ms4a6c" "Phf11b" "Ank"
## [271] "Nfib" "Kmo" "Myh4" "Lmnb1" "Tnfsf18" "Cox6a2"
## [277] "Epha5" "Kcnq1ot1" "Emp3" "Fabp4" "Cd7" "Kcnn2"
## [283] "Sh3tc2" "Gfod2" "Capg" "Cadm2" "Tk1" "Fpr2"
## [289] "Tsix" "Auts2" "Lockd" "Spc25" "Ly86" "Ntm"
## [295] "Klrb1b" "Amph" "Cxcl16" "Blk" "Ifit1bl1" "Timp2"DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB) +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB)
DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)
decide PCs to use
ElbowPlot(GEX.seur,ndims = 50)
PCs <- 1:20
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 20)## Computing nearest neighbor graph## Computing SNNGEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 1.2)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 25338
## Number of edges: 669263
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7603
## Number of communities: 20
## Elapsed time: 5 secondsRun UMAP/tSNE
GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 50, seed.use = 287)## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 11:40:13 UMAP embedding parameters a = 0.9922 b = 1.112## 11:40:13 Read 25338 rows and found 20 numeric columns## 11:40:13 Using Annoy for neighbor search, n_neighbors = 50## 11:40:13 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 11:40:16 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\Rtmpe8iaop\filee42c639d6904
## 11:40:16 Searching Annoy index using 1 thread, search_k = 5000
## 11:40:29 Annoy recall = 100%
## 11:40:30 Commencing smooth kNN distance calibration using 1 thread
## 11:40:33 Initializing from normalized Laplacian + noise
## 11:40:34 Commencing optimization for 200 epochs, with 1861048 positive edges
## 11:41:08 Optimization finishedDimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)
FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))
DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "FB.info", split.by = "FB.info",
ncol = 5, cols = color.FB)
DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "FB.info", split.by = "FB.info",
ncol = 5, cols = color.FB)
GEX.seur$cnt <- factor(gsub(".1$|.2$","",as.character(GEX.seur$FB.info)),
levels = c("P05.CTR","P05.MIG","P05.GFP",
"P28.CTR","P28.MIG","P28.GFP"))
color.cnt <- color.FB[c(10,9,7,
5,4,1)]DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
check some markers
markers.mig <- c("Ptprc",#"Cd3d","Cd3e","Cd19",
"Cd74","Lyz2","Ccl4",
"Aif1","P2ry12","C1qa","Spp1",
"Top2a","Pcna","Mki67","Mcm6",
"Cx3cr1","Il4ra","Il13ra1","Spp1-EGFP",
"Fabp5","Hmox1","Ms4a7","Cenpa")FeaturePlot(GEX.seur,
features = markers.mig,
ncol = 4)
DotPlot(GEX.seur, features = c("Cx3cr1","Spp1","Aif1", # Aif1: Iba1
"Fcer1g","Il4ra","Il13ra1"))
VlnPlot(GEX.seur, features = c("Cx3cr1","Spp1","Aif1", # Aif1: Iba1
"Fcer1g","Il4ra","Il13ra1"))
GEX.seur@meta.data$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
levels = c(1,2,7,5,
8,14,13,16,
9,
0,4,6,10,3,17,
11,12,
15,18,19))DotPlot(GEX.seur, features = markers.mig, group.by = "sort_clusters",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)
cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.info,
clusters=GEX.seur$sort_clusters)[c(3:12),],
main = "Cell Count",
gaps_row = c(2,4,5,7,9),
#gaps_col = c(12,13,14),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.info,
clusters=GEX.seur$sort_clusters)[c(3:12),]),
main = "Cell Ratio",
gaps_row = c(2,4,5,7,9),
#gaps_col = c(12,13,14),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
markers
# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "sort_clusters"
GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE,
min.pct = 0.05,
test.use = "MAST",
logfc.threshold = 0.25)
#GEX.markers.pre <- read.table("sc10x_LYN.marker.sort0811.csv", header = TRUE, sep = ",")
GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)## # A tibble: 160 x 7
## # Groups: cluster [20]
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
## 1 6.21e-297 0.422 0.995 0.988 1.14e-292 1 Rps11
## 2 1.47e-280 0.533 0.961 0.89 2.70e-276 1 Ltc4s
## 3 1.17e-199 0.453 0.966 0.928 2.15e-195 1 Maf
## 4 8.30e-180 0.548 0.816 0.677 1.52e-175 1 Crybb1
## 5 2.76e-174 0.529 0.817 0.722 5.07e-170 1 Ccr5
## 6 2.12e-146 0.502 0.654 0.528 3.89e-142 1 2410006H16Rik
## 7 1.84e-108 0.426 0.648 0.56 3.38e-104 1 Fchsd2
## 8 8.56e-105 0.426 0.428 0.33 1.57e-100 1 Fcgrt
## 9 0 0.812 0.91 0.669 0 2 Crybb1
## 10 1.30e-206 0.672 0.454 0.199 2.39e-202 2 Igfbp4
## # ... with 150 more rowsGEX.markers.sort <- GEX.markers.pre
GEX.markers.sort$cluster <- factor(as.character(GEX.markers.sort$cluster),
levels = levels(GEX.seur$sort_clusters))
markers.pre_t48 <- (GEX.markers.sort %>% group_by(cluster) %>%
filter(pct.1>0.1) %>%
top_n(n = 48, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene
markers.pre_t60 <- (GEX.markers.sort %>% group_by(cluster) %>%
filter(pct.1>0.1 & gene %in% grep("Rps|Rpl",GEX.markers.sort$gene,invert = T,value = T)) %>%
top_n(n = 60, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene
markers.pre_t120 <- (GEX.markers.sort %>% group_by(cluster) %>%
filter(pct.1>0.1 & gene %in% grep("Rps|Rpl",GEX.markers.sort$gene,invert = T,value = T)) %>%
top_n(n = 120, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$geneDotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[1:64])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[65:128])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[129:192])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[193:256])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[257:320])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[321:384])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[385:448])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[449:512])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[513:576])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[577:640])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[641:714])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
check.ref1 <- c("Tmem119","Hexb","Slc2a5","P2ry12",
"Siglech","Trem2", # microglia
"Mrc1","Lyve1","Cd163","Siglec1", # CAM
"Ly6c2","Ccr2","Plac8","Anxa8",
"Nr4a1", # Monocyte-derived
"Flt3","Zbtb46","Batf3","Itgae",
"Clec9a", #
"Tubb3","Actl6b" # locally addede neuron markers
)
DotPlot(GEX.seur, features = check.ref1, group.by = "sort_clusters",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)## Warning in FetchData(object = object, vars = features, cells = cells): The
## following requested variables were not found: Anxa8
C15: BAM (Brain Associated Macrophages)
C18: DC
C19: Neuron
others: Microglia
c("Xist","Tsix") %in% VariableFeatures(GEX.seur)## [1] TRUE TRUEc("Xist","Tsix","Rps4x") %in% VariableFeatures(GEX.seur)## [1] TRUE TRUE FALSEc("Rps4y1","Eif1ay","Ddx3y") %in% VariableFeatures(GEX.seur)## [1] FALSE FALSE FALSEc("Xist","Tsix","Uty","Ddx3y","Eif2s3y","Kdm5d") %in% VariableFeatures(GEX.seur)## [1] TRUE TRUE FALSE FALSE FALSE FALSEmarkers.Sci2019 <- c("Tmem119","Hexb","Slc2a5","P2ry12","Siglech","Trem2","Sparc","Olfml3", # Microglia
"Mrc1","Lyve1","Cd163","Siglec1","Pf4","Ms4a7","Cbr2", # CAMs
"Ly6c2","Ccr2","Plac8","Anxa8","Nr4a1","Mertk","Zbtb26","Cd209a", # Monocyte-derived
"Flt3","Zbtb46","Batf3","Itgae","Clec9a", # DCs
"Tubb3" # local added
)FeaturePlot(GEX.seur,
features = markers.Sci2019,
ncol = 4)## Warning in FetchData(object = object, vars = c(dims, "ident", features), : The
## following requested variables were not found: Anxa8
check.Cell2017 <- c("Hexb","Cst3","P2ry12","Cx3cr1",
"Cd9","Ctsd","Cst7","Lpl",
"Mrc1","Cd163","S100a4","Cd74",
"Cd79b","Rag1","Trbc2","Nkg7",
"S100a9","Camp"
)
DotPlot(GEX.seur, features = check.Cell2017, group.by = "sort_clusters",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)## Warning in FetchData(object = object, vars = features, cells = cells): The
## following requested variables were not found: Rag1, S100a9, Camp
DoubletFinder
library(DoubletFinder)## Loading required package: KernSmooth## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009## Loading required package: ROCR
## NULL## [1] "Creating 8446 artificial doublets..."
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."## Centering and scaling data matrix## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Computing pANN..."
## [1] "Classifying doublets.."## [1] "Creating 8446 artificial doublets..."
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Computing pANN..."
## [1] "Classifying doublets.."
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
ncol = 2,
pt.size = 0.05, group.by = "sort_clusters")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
ncol = 2,
pt.size = 0, group.by = "sort_clusters")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
ncol = 2,
pt.size = 0.05, group.by = "sort_clusters", split.by = "DoubletFinder0.05")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
ncol = 2,
pt.size = 0.05, group.by = "sort_clusters", split.by = "DoubletFinder0.1")
table(GEX.seur@meta.data[,c("sort_clusters","DoubletFinder0.05")])## DoubletFinder0.05
## sort_clusters Doublet Singlet
## 1 38 4348
## 2 286 3100
## 7 47 949
## 5 119 1613
## 8 65 825
## 14 15 236
## 13 17 385
## 16 12 151
## 9 71 645
## 0 122 4795
## 4 103 1883
## 6 93 1348
## 10 38 557
## 3 150 2004
## 17 12 115
## 11 30 448
## 12 33 382
## 15 12 164
## 18 2 110
## 19 2 13table(GEX.seur@meta.data[,c("sort_clusters","DoubletFinder0.1")])## DoubletFinder0.1
## sort_clusters Doublet Singlet
## 1 67 4319
## 2 450 2936
## 7 101 895
## 5 234 1498
## 8 137 753
## 14 50 201
## 13 94 308
## 16 39 124
## 9 111 605
## 0 204 4713
## 4 152 1834
## 6 143 1298
## 10 83 512
## 3 311 1843
## 17 63 64
## 11 41 437
## 12 75 340
## 15 98 78
## 18 77 35
## 19 4 11preAnno
# preAnno1
GEX.seur$preAnno1 <- as.character(GEX.seur$seurat_clusters)
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(1)] <- "MIG.a1"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(2)] <- "MIG.a2"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(7)] <- "MIG.a3"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(9)] <- "MIG.a4"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(5)] <- "MIG.a5"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(8)] <- "MIG.a6"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(14)] <- "MIG.a7"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(13)] <- "MIG.a8"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(16)] <- "MIG.a9"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(0)] <- "MIG.b1"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(11)] <- "MIG.b2"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(4)] <- "MIG.b3"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(6)] <- "MIG.b4"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(12)] <- "MIG.b5"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(10)] <- "MIG.b6"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(3)] <- "MIG.b7"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(17)] <- "MIG.b8"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(18)] <- "DC"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(15)] <- "BAM"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(19)] <- "Neuron"
GEX.seur$preAnno1 <- factor(GEX.seur$preAnno1,
levels = c(paste0("MIG.a",1:9),
paste0("MIG.b",1:8),
"DC","BAM","Neuron"))
# preAnno2
GEX.seur$preAnno2 <- as.character(GEX.seur$preAnno1)
GEX.seur$preAnno2 <- gsub("1|2|3|4|5|6|7|8|9","",GEX.seur$preAnno2)
GEX.seur$preAnno2 <- factor(GEX.seur$preAnno2,
levels = c("MIG.a","MIG.b",
"DC","BAM","Neuron"))
# preAnno
GEX.seur$preAnno <- as.character(GEX.seur$preAnno2)
GEX.seur$preAnno <- gsub(".a$|.b$","",GEX.seur$preAnno)
GEX.seur$preAnno <- factor(GEX.seur$preAnno,
levels = c("MIG","DC","BAM","Neuron"))pp.raw.1a <- DimPlot(GEX.seur, reduction = "umap", label = T,label.size = 2.4, group.by = "preAnno1") +
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "sort_clusters")
pp.raw.1a
pp.raw.1b <- DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
pp.raw.1b
pp.raw.1c <- DimPlot(GEX.seur, reduction = "umap", group.by = "cnt", split.by = "cnt", cols = color.cnt, ncol = 3)
pp.raw.1c 
DimPlot(GEX.seur, reduction = "umap", label = T,label.size = 2.4, group.by = "preAnno1") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno2")
#saveRDS(GEX.seur, "GEX0811.seur.preAnno.rds")cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Spp1", reduction = "pca",dims = 1:2),
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB),
rel_widths = c(5,5.6),
ncol = 2)
pp.test <- FeaturePlot(GEX.seur, features = "Spp1", reduction = "pca",dims = 1:2, split.by = "FB.info", combine = F)
cowplot::plot_grid(
plotlist = pp.test,
ncol = 5
)
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", split.by = "FB.info", cols = color.FB, ncol = 5)
cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Spp1", reduction = "pca",dims = 3:4),
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB),
rel_widths = c(5,5.6),
ncol = 2)
cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$cnt,
clusters=GEX.seur$preAnno),
main = "Cell Count",
#gaps_row = c(2,4,5,7,9),
#gaps_col = c(12,13,14),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$cnt,
clusters=GEX.seur$preAnno)),
main = "Cell Ratio",
#gaps_row = c(2,4,5,7,9),
#gaps_col = c(12,13,14),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
table(GEX.seur@meta.data[,c("cnt","preAnno")])## preAnno
## cnt MIG DC BAM Neuron
## P05.CTR 2598 2 12 5
## P05.MIG 5362 11 32 5
## P05.GFP 5126 1 32 2
## P28.CTR 2372 3 14 0
## P28.MIG 4980 36 27 2
## P28.GFP 4597 59 59 1rowRatio(table(GEX.seur@meta.data[,c("cnt","preAnno")]))## preAnno
## cnt MIG DC BAM Neuron
## P05.CTR 0.9927397784 0.0007642339 0.0045854031 0.0019105846
## P05.MIG 0.9911275416 0.0020332717 0.0059149723 0.0009242144
## P05.GFP 0.9932183685 0.0001937609 0.0062003488 0.0003875218
## P28.CTR 0.9928840519 0.0012557555 0.0058601925 0.0000000000
## P28.MIG 0.9871159564 0.0071357780 0.0053518335 0.0003964321
## P28.GFP 0.9747667515 0.0125106022 0.0125106022 0.0002120441marker.fig0921 <- c("Ptprc","P2ry12","Spp1","Spp1-EGFP",
"Cx3cr1","Csf1r","Aif1","C1qa",
"Ccl4","Lpl","Fabp5","Csf1",
"Ctsf","Ccr5","Arsb","Itgb1",
"Cd81","Lag3","Cd34","H2-K1",
"Top2a","Pcna","Mki67","Mcm6",
"Cd74","H2-Aa","Ccr2","Cd209a",
"Mrc1","Pf4","Dab2","Lyve1",
"Tubb3","Actl6b","Syt1","Sst"
)
pp.raw.2a <- FeaturePlot(GEX.seur, features = marker.fig0921, ncol = 4,raster = T, pt.size = 2.75)
pp.raw.2a
pp.raw.2b <- DotPlot(GEX.seur, features = marker.fig0921, group.by = "preAnno1",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)
pp.raw.2b